洛谷P4491 [HAOI2018]染色

设为钦定种颜色，每种颜色恰好出现了次的方案数，设为恰好种颜色，每种颜色恰好出现了次的方案数

考虑钦定k种颜色能得到的方案数，可得

阶乘预处理后即可时间内算出

然后考虑和之间的关系，有

二项式反演后可得

差卷积后乘上对应值即可得出答案，总复杂度

#include <bits/stdc++.h>
#include <ext/pb_ds/tree_policy.hpp>
#include <ext/pb_ds/trie_policy.hpp>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/priority_queue.hpp>

using namespace std;
using namespace __gnu_pbds;
using namespace __gnu_cxx;
typedef long long ll;
//tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> tr;
//__gnu_pbds::priority_queue<int, greater<int>, pairing_heap_tag> qu;
//typedef trie<string, null_type, trie_string_access_traits<>, pat_trie_tag, trie_prefix_search_node_update> pref_trie;
const int N = 10000005;
int n, m, k;
namespace polybase {//范围为1e9需要先取模,别忘记改模数
    const ll mod = 1004535809;
    int limit;
    ll _wn[25];

    ll fpow(ll x, ll r)
    {
        ll result = 1;
        while (r)
        {
            if (r & 1)result = result * x % mod;
            r >>= 1;
            x = x * x % mod;
        }
        return result;
    }

    int _ = []
    {
        for (int i = 0; i <= 23; i++)_wn[i] = fpow(3, (mod - 1) >> i);
        return 0;
    }();

    inline int norm(int n) { return 1 << __lg(n * 2 - 1); }

    void NTT(ll *A, int type)
    {
        int i, j = limit >> 1, k, l, c = 0;
        ll u, v, w, wn;
        for (i = 1; i < limit - 1; i++)
        {
            if (i < j)swap(A[i], A[j]);
            for (k = limit >> 1; (j ^= k) < k; k >>= 1);
        }

        for (l = 2; l <= limit; l <<= 1)
        {
            i = l >> 1, wn = _wn[++c];
            for (j = 0; j < limit; j += l)
            {
                w = 1;
                for (k = j; k < j + i; k++)
                {
                    u = A[k], v = A[k + i] * w % mod;
                    A[k] = u + v;
                    if (A[k] >= mod)A[k] -= mod;

                    A[k + i] = u - v;
                    if (A[k + i] < 0)A[k + i] += mod;
                    w = w * wn % mod;
                }
            }
        }
        if (type == -1)
        {
            ll inv = fpow(limit, mod - 2);
            for (i = 0; i < limit; i++)A[i] = A[i] * inv % mod;
            for (i = 1; i < limit / 2; i++)swap(A[i], A[limit - i]);
        }
    }

    struct poly : public vector<ll>
    {
        using vector<ll>::vector;
#define T (*this)

        poly modxk(int k) const
        {
            k = min(k, (int) size());
            return poly(begin(), begin() + k);
        }

        poly rev() const { return poly(rbegin(), rend()); }

        friend void NTT(poly &a, const int type) { NTT(a.data(), type); }

        friend poly operator*(const poly &x, const poly &y)
        {
            if (x.empty() || y.empty())return poly();
            poly a(x), b(y);
            int len = a.size() + b.size() - 1;
            limit = norm(len);
            a.resize(limit), b.resize(limit);
            NTT(a, 1);
            NTT(b, 1);
            for (int i = 0; i < limit; i++)
                a[i] = a[i] * b[i] % mod;
            NTT(a, -1);
            a.resize(len);
            return a;
        }

        poly operator+(const poly &b)
        {
            poly a(T);
            if (a.size() < b.size())
                a.resize(b.size());
            for (int i = 0; i < b.size(); i++)//用b.size()防止越界
            {
                a[i] += b[i];
                if (a[i] >= mod)a[i] -= mod;
            }
            return a;
        }

        poly operator-(const poly &b)
        {
            poly a(T);
            if (a.size() < b.size())
                a.resize(b.size());
            for (int i = 0; i < b.size(); i++)
            {
                a[i] -= b[i];
                if (a[i] < 0)a[i] += mod;
            }
            return a;
        }

        poly operator*(const ll p)
        {
            poly a(T);
            for (auto &x:a)
                x = x * p % mod;
            return a;
        }

        poly &operator<<=(int r) { return insert(begin(), r, 0), T; }//注意逗号,F(x)*(x^r)

        poly operator<<(int r) const { return poly(T) <<= r; }

        poly operator>>(int r) const { return r >= size() ? poly() : poly(begin() + r, end()); }

        poly &operator>>=(int r) { return T = T >> r; }//F[x]/(x^r)

        poly deriv()//求导
        {
            if (empty())return T;
            poly a(size() - 1);
            for (int i = 1; i < size(); i++)//注意是size()
                a[i - 1] = T[i] * i % mod;
            return a;
        }

        poly integ()//积分
        {
            poly a(size() + 1);
            for (int i = 1; i < a.size(); i++)//注意是a.size()
                a[i] = T[i - 1] * fpow(i, mod - 2) % mod;
            return a;
        }

        poly inv(int n)
        {
            poly a{fpow(T[0], mod - 2)};
            int k = 1;
            while (k < n)
            {
                k <<= 1;
                a = (a * 2 - modxk(k) * a * a).modxk(k);
            }
            return a.modxk(n);
        }

        poly sqrt(int n)//f[0]必须等于1
        {
            poly a{1};
            int k = 1;
            const ll inv2 = fpow(2, mod - 2);
            while (k < n)
            {
                k <<= 1;
                a = ((modxk(k) * a.inv(k)).modxk(k) + a) * inv2;
            }
            return a.modxk(n);
        }

        poly ln(int n)//需要保证f[0]=1
        {
            return (deriv() * inv(n)).integ().modxk(n);
        }

#undef T
    };
}
using namespace polybase;
ll fac[N], ifac[N];

int main()
{
    int p, q, u, v, w, x, y, z, T;
    fac[0] = 1;
    for (int i = 1; i <= N - 5; i++)
        fac[i] = fac[i - 1] * i % mod;
    ifac[N - 5] = fpow(fac[N - 5], mod - 2);
    for (int i = N - 5; i; i--)
        ifac[i - 1] = ifac[i] * i % mod;
    int s;
    cin >> n >> m >> s;
    poly W(m + 1);
    for (int i = 0; i <= m; i++)
        scanf("%d", &W[i]);
    poly F(m + 1);
    for (int i = min(n / s, m); i >= 0; i--)
        F[i] = fac[m] * ifac[m - i] % mod * fac[n] % mod * fpow(ifac[s], i) % mod *
               ifac[n - i * s] % mod * fpow(m - i, n - i * s) % mod;
    poly H(m + 1);
    for (int i = 0; i <= m; i++)
        H[m - i] = (i & 1 ? mod - 1 : 1) * ifac[i] % mod;
    poly G = F * H;
    ll ans = 0;
    for (int i = 0; i <= m; i++)
        ans = (ans + G[m + i] * ifac[i] % mod * W[i]) % mod;
    cout << ans;
    return 0;
}